	********************************************************************************
	****
	**** Date: 04.25.2022							 
	**** Update: 07.12.2023					 
	**** Author: JW josephgwright@gmail.com 				 	
	**** NOTE: This program has been executed in Stata 18.0	 
	**** 
	**** Using data files in $dir: 
	****
	****  	pers-use.dta
	**** 	eap.dta
	****
	*********************************************************************************
	
			
		capture log close
		log using Ch4.log, replace

	******************************
	**** Set directory, seed *****
	******************************
		set more off 
		set matsize 1000
		global seed ="984353"
		set scheme plotplain
		cd "$dir"

		* Load data *
		use pers-use,clear
		drop if persparty==.
		
		* GNB 2012: non-personalist leader dies a naturual death 9 Jan 2012; 
		* interim president takes over and their is a coup in April 
		* drop the interim leader who is unelected *
 		drop if GNB==1
		
		* Erosion indicators *
		forval i = 5(1)15 {
 			di `i'
			local c = `i'/100
 			gen decline`i' = (v2x_polyarchy - ivdem)<=-`c'  if  ivdem~=. & v2x_polyarchy~=.
			xtset lid year 
			gen ld`i'=l.decline`i'  
		}
		
		* Update for collapse in 2020 *
		recode gwf_back (0=1) if year==2020 & (country=="Mali")
				
		* Generate some confounders *
		gen econcrisis = l12gr<-2 if l12gr~=.
		tab econcrisis	
		gen d1 =gwf_duration
		gen d2 = d1^2
		gen d3 = d1^3
		replace lnparty = ln(1+partyage)
		gen osupdem = l1supdem if year==min
		replace osupdem = l2supdem if year==min & osupdem==.
		replace osupdem = supdem if year==min & osupdem==.
		egen isupdem = max(osupdem),by(lid)
		gen opolar = l1polar  if year==min
		replace opolar  = l2polar  if year==min & opolar ==.
		replace opolar  = polar  if year==min & opolar ==.
		egen ipolar  = max(opolar),by(lid)
		gen time= year-1990
		
		* Missingness in support for democracy *
		gen miss = isupdem==.
		local var = "gwf_back ld ivdem lpop year persparty create"
		foreach v of local var {
			ttest `v',by(miss)unequal
		} 
 	 
		* Sample *
		qui reg v2x_polyarchy ld ivdem
		sum v2x_polyarchy $ldv year if e(sample)==1
		sum v2x_polyarchy $ldv year  
		
		
		*********************
		** Civil liberites **
		*********************
		qui reg v2x_polyarchy  create ld  if persparty~=.
			gen sample = e(sample)==1
			egen c=count(year) if e(sample)==1,by(lid) 	
			hist c
			gen election = v2xel_elecparl==1 | v2xel_elecpres==1
			tsset cowcode year
			gen f1election  =f.election
			gen l1election  =l.election
			
		local var = "v2x_clpol v2x_clpriv v2x_clphy frepress"
		foreach v of local var {
			 qui gen o`v'=(l1`v') if minyr==year
			 qui egen i`v'=max(o`v'),by(lid)
		}
		* flip scales to measure more repression and  * v2x_clphy already flipped; v2caviol measures more violence
		local var = "v2x_clpol v2x_clpriv"
		foreach v of local var {
			replace `v'=`v'*-1
		}
		* rescale rescale on 0,1 *
		local var = "v2x_clpol v2x_clpriv v2x_clphy frepress"
		foreach v of local var {
			qui sum `v'
			replace `v'=`v'+abs(r(min))
			qui sum `v'
			replace `v'=`v'/r(max)
		}
		sum v2x_clpol v2x_clpriv v2x_clphy v2caviol
		corr v2x_clpol v2x_clpriv v2x_clphy v2caviol
		swilk v2x_clpol v2x_clpriv v2x_clphy v2caviol
		
		local var = "v2x_clpol v2x_clpriv"
		foreach v of local var {
			replace `v'=`v'^(1/3)
			qui sum `v'
			replace `v'=(`v'-r(mean))/r(sd)
		}	
		local var = "v2x_clphy"
		foreach v of local var {
			replace `v'=`v'^(1/12)
			qui sum `v'
			replace `v'=(`v'-r(mean))/r(sd)
		}	
		local var = "frepress"
		foreach v of local var {
		qui sum `v'
			replace `v'=(`v'-r(mean))/r(sd)
		}
		swilk v2x_clpol v2x_clpriv v2x_clphy frepress
		sum v2x_clpol v2x_clpriv v2x_clphy frepress
		


		  * Differenced Lag DV *
		  reghdfe d.v2x_clpol persparty ld ivdem l1v2x_clpol l2v2x_clpol time,a(cowcode)cluster(lid) 
			    qui xtreg d.v2x_clpol time persparty ld ivdem l1v2x_clpol l2v2x_clpol,i(cowcode)fe cluster(cowcode) 
			    qui predict e_residuals_1, e
				xtistest e_res, lags(2)
				drop e_res
		  * FE with initial level *
		   reghdfe v2x_clpol persparty ld ivdem iv2x_clpol time,a(cowcode)cluster(lid)
			    qui xtreg v2x_clpol time persparty ld ivdem iv2x_clpol,cluster(cowcode)fe
			    qui predict e_residuals_1, e
				xtistest e_res, lags(2)
				drop e_res
		  * Rescale to get estimates that fit scale of treatment -- only for visualizing coefficient estimates 
			replace ivdem = ivdem*10
			replace ipi = ipi*10
			replace l1v2x_clpol=l1v2x_clpol*50
			replace l2v2x_clpol=l2v2x_clpol*50
		
		  * Dynamic panel baseline model *
		   reghdfe v2x_clpol persparty ld ivdem l1v2x_clpol l2v2x_clpol,a(cowcode)cluster(lid)
		   est store civlib1
				qui xtreg v2x_clpol time persparty ld ivdem l1v2x_clpol l2v2x_clpol,cluster(cowcode)fe
				qui predict e_residuals_1, e
				xtistest e_res, lags(2)
				drop e_res	
				
			gen newparty  =(partyage^(1/2))*-1
			qui sum newparty
			replace newparty = newparty +abs(r(min))
			qui sum newparty 
			replace newparty = newparty/r(max)
			sum newparty
			hist newparty
			swilk newparty lnparty partyage
			
			reghdfe v2x_clpol persparty time ld ivdem newparty l1v2x_clpol l2v2x_clpol,a(cowcode)cluster(lid)
			est store civlib2
			reghdfe v2x_clpol persparty time ld ivdem ipi l1v2x_clpol l2v2x_clpol,a(cowcode)cluster(lid)
			est store civlib3
			reghdfe v2x_clpol persparty time ld ivdem i_popul l1v2x_clpol l2v2x_clpol,a(cowcode)cluster(lid)
			est store civlib4
		    reghdfe v2x_clpol persparty time ld ivdem l1v2x_clpol l2v2x_clpol if year<2020,a(cowcode)cluster(lid)

			label var persparty "{bf:Party personalism}"
			label var ld "Democracy age"
			label var ivdem `""Initial democracy" "level          ""'
			label var newparty "New parties index"
			label var ipi `""Initial party     " "institutionalization""'
			label var i_popul "Party populism"
			coefplot (civlib1, msymbol(d))(civlib2, msymbol(P)) (civlib3, msymbol(T)) (civlib4, msymbol(Oh)), order(persparty)  ///
				drop(_cons l1v2x_clpol l2v2x_clpol time) xline(0) msymbol(d) mfcolor(white) grid(glcolor(gs15)) ///
				levels(95 90) legend(lab(3 "Baseline")lab(6 "+ New party index") ///
				lab(9 "+ Initial party institutionalization")lab(12 "+ Party populism")  order(3 6 9 12) ///
				size(small) pos(6) col(2) ring(1)) xsize(2) ysize(2) xlab(-.2(.1).3)  ///
				xtitle("        Coefficient estimate", size(small))  ///
				ciopts(lwidth(thin)) aspectratio(1.1) scale(.75) title(Government repression of civil liberties, size(medium) height(2))
			gr export "$dir\golden\Ch4-PersParty-Civliberties.pdf",as(pdf)replace 
		  
		  * Create party * 
		   reghdfe v2x_clpol create ld ivdem l1v2x_clpol l2v2x_clpol,a(cowcode year)cluster(lid)
		   est store c1
		   
		   * IFE *
			regife v2x_clpol persparty ld ivdem l1v2x_clpol l2v2x_clpol,a(cowcode year)ife(cowcode year,1)vce(cluster lid)
			est store ife1
			regife v2x_clpol persparty ld ivdem l1v2x_clpol l2v2x_clpol,a(cowcode year)ife(cowcode year,2)vce(cluster lid)	
			est store ife2
			label var l1v2x_clpol  `""Civil     " "liberties{sub:t-1}""'
			label var l2v2x_clpol `""Civil     " "liberties{sub:t-2}""'
			label var create   `""Create" "Party ""'
			coefplot (civlib1, msymbol(d))(c1, msymbol(P)) (ife1, msymbol(T)) (ife2, msymbol(Oh)), order(persparty create)  ///
				drop(_cons  time) xline(0) msymbol(d) mfcolor(white) grid(glcolor(gs15)) ///
				levels(95 90) legend(lab(3 "Party personalism, dynamic panel")lab(6 "Create party, dynamic panel") ///
				lab(9 "IFE, 1")lab(12 "IFE, 2")  order(3 6 9 12) ///
				size(small) pos(6) col(2) ring(1)) xsize(2) ysize(2) xlab(-.2(.1).3)  ///
				xtitle("        Coefficient estimate", size(small))  ///
				ciopts(lwidth(thin)) aspectratio(1.1) scale(.75) title(Government repression of civil liberties, size(medium) height(2))
			gr export "$dir\golden\T-Civlib-Robust.pdf",as(pdf)replace 
				
		   * Check dynamic panel estimate against confounders *
		    global d = "persparty"
			gen beta=.
			gen hi90=.
			gen lo90=.
			gen n =_n
			gen hi=.
			gen lo=.
			gen varname = " "
		   	local i =1
			local var = "l1v2xps_party ipi i_popu election i.v2elparlel v2pavote v2paseats ipolar isupdem priormil pres econcrisis gdp lpop oilgas"
			foreach v of local var {
				di "`v'"
				xtset lid year
				reghdfe v2x_clpol `v' persparty time ld ivdem l1v2x_clpol l2v2x_clpol,a(cowcode)cluster(lid) 
				lincom $d
				qui nlcom _b[$d],post
				matrix beta =e(b)  
				local b = beta[1,1]
				qui replace beta=`b' if n==`i'
				matrix var = e(V) 
				local se =var[1,1]
				qui replace hi = `b' + sqrt(`se')*1.96 if n==`i'
				qui replace lo = `b' - sqrt(`se')*1.96 if n==`i'
				qui replace hi90 = `b' + sqrt(`se')*1.645 if n==`i'
				qui replace lo90 = `b' - sqrt(`se')*1.645 if n==`i'
				qui replace varname = "`v'" if n==`i'
				local i = `i' +1
			 }
			label define varlab  1 "Party instit." 2 "Initial party inst." 3 "Populism" ///
				4 "Election" 5 "Electoral system" 6 `""Ruling party" "leg. seat share""' ///
				7 `""Ruling party" "vote share""' 8 "Polarization" 9 `""Citizen support" "for democracy""' 10 "Prior military" ///
				11 "Presidential" 12 "Economic crisis" 13 "GDP per capita" 14 "Population" 15 "Oil rents",replace
			label values n varlab
			twoway (scatter beta n if n<=15,mcol(blue)yscale(range(-.01 0.08))yline(.1170501,lcol(gs4)lpat(dash_dot))) ///
				(rspike hi lo n if n<=15,lw(vthin)lcol(blue)ylab(0(.1)0.2)) ///
				(rspike hi90 lo90 n if n<=15,lcol(blue)lw(medium)ytitle("{&beta}{sub:Party personalism}", ///
				size(large)height(4))tit(Civil liberties repression)subtit("Covariate adjustment, dynamic panel model",size(small))   ///
				xtitle(Added covariate,height(-10))yline(0,lpat(dash)lcol(red))xlab(1(1)15,valuelabel angle(90))legend(off) ///
				note("Baseline covariates: Democracy age, initial level of democracy and two lags of the outcome variable",  ///
				size(vsmall)pos(6)))
			gr export "$dir\golden\T-Civlib-Covariates.pdf",as(pdf)replace 
			
			qui reghdfe v2x_clpol persparty time ld ivdem l1v2x_clpol l2v2x_clpol,a(cowcode)cluster(lid) 
			lincom persparty
			sum persparty isupdem if e(sample)==1
			qui reghdfe v2x_clpol persparty time ld ivdem l1v2x_clpol l2v2x_clpol if isupdem~=.,a(cowcode)cluster(lid)
			lincom persparty
			qui reghdfe v2x_clpol isupdem persparty time ld ivdem l1v2x_clpol l2v2x_clpol if isupdem~=.,a(cowcode)cluster(lid) 
			lincom persparty

 			gen misssupdem = isupdem==.
			ttest gwf_back, by(misssupdem)
 			 
		*******************************		
		** T-tests by party creation **
		*******************************		
			gen e=.	
			ttest v2x_poly,by(create)
			local  m1=r(mu_1) 
			local se1 = r(sd_1)/(sqrt(r(N_1)))
			local  m2=r(mu_2)
			local se2 = r(sd_2)/(sqrt(r(N_2)))
			replace e=`m1' if _n==1
			replace hi = `m1' + 1.96* `se1' if _n==1
			replace lo = `m1' - 1.96* `se1'  if _n==1
			replace e=`m2' if _n==2
			replace hi = `m2' + 1.96*`se2' if _n==2
			replace lo = `m2' - 1.96*`se2' if _n==2				
			twoway (bar e n if n<=2,barwidth(.5)bcol(gs13)ytit(Democracy level)saving(h1.gph,replace)) ///
				(rspike hi lo n if n<=2,ylab(.5(.1).8)col(gs1)legend(off)xtit("")tit(Democracy decay) ///
				xlab(1 "No party creation" 2 "Leader creates party")xscale(range(.8 2.2)))
			ttest decline10,by(create)
			local  m1=r(mu_1) 
			local se1 = r(sd_1)/(sqrt(r(N_1)))
			local  m2=r(mu_2)
			local se2 = r(sd_2)/(sqrt(r(N_2)))
			replace e=`m1' if _n==1
			replace hi = `m1' + 1.96* `se1' if _n==1
			replace lo = `m1' - 1.96* `se1'  if _n==1
			replace e=`m2' if _n==2
			replace hi = `m2' + 1.96*`se2' if _n==2
			replace lo = `m2' - 1.96*`se2' if _n==2				
			twoway (bar e n if n<=2,barwidth(.5)bcol(gs13)ytit("Probability of democratic decline of 10 % or more")saving(h2.gph,replace)) ///
				(rspike hi lo n if n<=2,ylab(0(.02).1)col(gs1)legend(off)xtit("")tit(Democratic erosion) ///
				xlab(1 "No party creation" 2 "Leader creates party")xscale(range(.8 2.2)))
			ttest gwf_back,by(create)
			local  m1=r(mu_1) 
			local se1 = r(sd_1)/(sqrt(r(N_1)))
			local  m2=r(mu_2)
			local se2 = r(sd_2)/(sqrt(r(N_2)))
			replace e=`m1' if _n==1
			replace hi = `m1' + 1.95* `se1' if _n==1
			replace lo = `m1' - 1.95* `se1'  if _n==1
			replace e=`m2' if _n==2
			replace hi = `m2' + 1.95*`se2' if _n==2
			replace lo = `m2' - 1.95*`se2' if _n==2				
			twoway (bar e n if n<=2,barwidth(.5)bcol(gs13)ytit(Probability of democratic collapse)saving(h3.gph,replace)) ///
				(rspike hi lo n if n<=2,ylab(0(.02).1)col(gs1)legend(off)xtit("")tit(Democratic collapse) ///
				xlab(1 "No party creation" 2 "Leader creates party")xscale(range(.8 2.2)))
			gr combine h1.gph h2.gph h3.gph, xsize(4)ysize(2)col(3)
			gr export "$dir\golden\Ch4-Democracy-ttests.pdf",as(pdf)replace 
	 
			*** Initial support for democracy ***
			qui centile isupdem,centile(50)
			local c = r(c_1)
			tab gwf_back create if isupdem<`c'  & isupdem~=.,col
			ttest gwf_back if isupdem<`c'  & isupdem~=.,by(create)level(90)
			local  m1=r(mu_1) 
			local se1 = r(sd_1)/(sqrt(r(N_1)))
			local  m2=r(mu_2)
			local se2 = r(sd_2)/(sqrt(r(N_2)))
			replace e=`m1' if _n==1
			replace hi = `m1' + 1.95* `se1' if _n==1
			replace lo = `m1' - 1.95* `se1'  if _n==1
			replace e=`m2' if _n==2
			replace hi = `m2' + 1.95*`se2' if _n==2
			replace lo = `m2' - 1.95*`se2' if _n==2				
			twoway (bar e n if n<=2,barwidth(.5)bcol(gs13)ytit(Probability of democratic collapse)saving(h1.gph,replace)) ///
				(rspike hi lo n if n<=2,ylab(0(.02).06)yscale(range(0 .07)) ///
				col(gs1)legend(off)xtit("")tit({bf:Low} initial support for democracy) ///
				xlab(1 "No create party (5)" 2 "Create party (7)")xscale(range(.8 2.2)))
			tab gwf_back create  if isupdem>`c'  & isupdem~=.,col
			ttest gwf_back if isupdem>`c' &  isupdem~=.,by(create)level(90)
			local  m1=r(mu_1) 
			local se1 = r(sd_1)/(sqrt(r(N_1)))
			local  m2=r(mu_2)
			local se2 = r(sd_2)/(sqrt(r(N_2)))
			replace e=`m1' if _n==1
			replace hi = `m1' + 1.95* `se1' if _n==1
			replace lo = `m1' - 1.95* `se1'  if _n==1
			replace e=`m2' if _n==2
			replace hi = `m2' + 1.95*`se2' if _n==2
			replace lo = `m2' - 1.95*`se2' if _n==2				
			twoway (bar e n if n<=2,barwidth(.5)bcol(gs13)ytit(Probability of democratic collapse)saving(h2.gph,replace)) ///
				(rspike hi lo n if n<=2,ylab(0(.02).06)yscale(range(0 .07)) ///
				col(gs1)legend(off)xtit("")tit({bf:High} initial support for democracy) ///
				xlab(1 "No create party (0)" 2 "Create party (3)")xscale(range(.8 2.2)))
			tab gwf_back create  if  isupdem==.,col
			ttest gwf_back if isupdem==.,by(create)level(90)
			local  m1=r(mu_1) 
			local se1 = r(sd_1)/(sqrt(r(N_1)))
			local  m2=r(mu_2)
			local se2 = r(sd_2)/(sqrt(r(N_2)))
			replace e=`m1' if _n==1
			replace hi = `m1' + 1.95* `se1' if _n==1
			replace lo = `m1' - 1.95* `se1'  if _n==1
			replace e=`m2' if _n==2
			replace hi = `m2' + 1.95*`se2' if _n==2
			replace lo = `m2' - 1.95*`se2' if _n==2				
			twoway (bar e n if n<=2,barwidth(.5)bcol(gs13)ytit(Probability of democratic collapse)saving(h3.gph,replace)) ///
				(rspike hi lo n if n<=2,ylab(0(.02).06)yscale(range(0 .07)) ///
				col(gs1)legend(off)xtit("")tit({bf:Missing data} on support for democracy) ///
				xlab(1 "No create party (9)" 2 "Create party (11)")xscale(range(.8 2.2)))
			gr combine h1.gph h2.gph h3.gph,xsize(8)tit(Democratic collapse)col(3)
			gr export "$dir\golden\T-Democracy-ttests-by-demsupport.pdf",as(pdf)replace 
			
			*** Party populism **
			qui centile i_pop,centile(50)
			local c = r(c_1)
			tab gwf_back create if i_pop<`c'  & i_pop~=.,col
			ttest gwf_back if i_pop<`c'  & i_pop~=.,by(create)level(90)
			local  m1=r(mu_1) 
			local se1 = r(sd_1)/(sqrt(r(N_1)))
			local  m2=r(mu_2)
			local se2 = r(sd_2)/(sqrt(r(N_2)))
			replace e=`m1' if _n==1
			replace hi = `m1' + 1.645* `se1' if _n==1
			replace lo = `m1' - 1.645* `se1'  if _n==1
			replace e=`m2' if _n==2
			replace hi = `m2' + 1.645*`se2' if _n==2
			replace lo = `m2' - 1.645*`se2' if _n==2
			twoway (bar e n if n<=2,barwidth(.5)bcol(gs13)ytit(Probability of democratic collapse)saving(h1.gph,replace)) ///
				(rspike hi lo n if n<=2,ylab(0(.02).06)yscale(range(0 .06)) ///
				col(gs1)legend(off)xtit("")tit("{bf:Low} populism")text(.05 1.5 "p<0.01") ///
				xlab(1 "Not create party (2)" 2 "Create party (6)")xscale(range(.8 2.2)))
			tab gwf_back create if i_pop>=`c'  & i_pop~=.,col
			ttest gwf_back if i_pop>=`c'  & i_pop~=.,by(create)level(90)
			local  m1=r(mu_1) 
			local se1 = r(sd_1)/(sqrt(r(N_1)))
			local  m2=r(mu_2)
			local se2 = r(sd_2)/(sqrt(r(N_2)))
			replace e=`m1' if _n==1
			replace hi = `m1' + 1.645* `se1' if _n==1
			replace lo = `m1' - 1.645* `se1'  if _n==1
			replace e=`m2' if _n==2
			replace hi = `m2' + 1.645*`se2' if _n==2
			replace lo = `m2' - 1.645*`se2' if _n==2				
			twoway (bar e n if n<=2,barwidth(.5)bcol(gs13)ytit(Probability of democratic collapse)saving(h2.gph,replace)) ///
				(rspike hi lo n if n<=2,ylab(0(.02).06)yscale(range(0 .06)) ///
				col(gs1)legend(off)xtit("")tit("{bf:High} populism")text(.05 1.5 "p<0.05") ///
				xlab(1 "Not create party (11)" 2 "Create party (13)")xscale(range(.8 2.2)))
			
			*** Party age ***
			egen byear = min(year),by(lid)
			replace minyr = byear
			gen oage = partyage if year==minyr
			egen iage=max(oage),by(lid)
			centile iage if year==min,centile(35)
			local c = 10
			tab gwf_back create if iage<=`c' ,col
			ttest gwf_back if iage<=`c',by(create)level(90)
			local  m1=r(mu_1) 
			local se1 = r(sd_1)/(sqrt(r(N_1)))
			local  m2=r(mu_2)
			local se2 = r(sd_2)/(sqrt(r(N_2)))
			replace e=`m1' if _n==1
			replace hi = `m1' + 1.645* `se1' if _n==1
			replace lo = `m1' - 1.645* `se1'  if _n==1
			replace e=`m2' if _n==2
			replace hi = `m2' + 1.645*`se2' if _n==2
			replace lo = `m2' - 1.645*`se2' if _n==2				
			twoway (bar e n if n<=2,barwidth(.5)bcol(gs13)ytit(Probability of democratic collapse)saving(h3.gph,replace)) ///
				(rspike hi lo n if n<=2,ylab(0(.02).06)yscale(range(0 .06)) ///
				col(gs1)legend(off)xtit("")tit({bf:Young} parties)text(.05 1.5 "p<0.10") ///
				xlab(1 "No create party (6)" 2 "Create party (17)")xscale(range(.8 2.2)))
			tab gwf_back create  if iage>`c'  ,col
			ttest gwf_back if iage>`c'  ,by(create)level(90)
			local  m1=r(mu_1) 
			local se1 = r(sd_1)/(sqrt(r(N_1)))
			local  m2=r(mu_2)
			local se2 = r(sd_2)/(sqrt(r(N_2)))
			replace e=`m1' if _n==1
			replace hi = `m1' + 1.645* `se1' if _n==1
			replace lo = `m1' - 1.645* `se1'  if _n==1
			replace e=`m2' if _n==2
			replace hi = `m2' + 1.645*`se2' if _n==2
			replace lo = `m2' - 1.645*`se2' if _n==2				
			twoway (bar e n if n<=2,barwidth(.5)bcol(gs13)ytit(Probability of democratic collapse)saving(h4.gph,replace)) ///
				(rspike hi lo n if n<=2,ylab(0(.02).06)yscale(range(0 .06)) ///
				col(gs1)legend(off)xtit("")tit({bf:Old} parties)text(.05 1.5 "p<0.05") ///
				xlab(1 "No create party (8)" 2 "Create party (4)")xscale(range(.8 2.2)))
 			
			*** Party system institutionalization ***
			centile ipi if year==min,centile(30)
			local c = r(c_1)
			tab gwf_back create if ipi<`c' ,col
			ttest gwf_back if ipi<`c',by(create)level(95)
			local  m1=r(mu_1) 
			local se1 = r(sd_1)/(sqrt(r(N_1)))
			local  m2=r(mu_2)
			local se2 = r(sd_2)/(sqrt(r(N_2)))
			replace e=`m1' if _n==1
			replace hi = `m1' + 1.645* `se1' if _n==1
			replace lo = `m1' - 1.645* `se1'  if _n==1
			replace e=`m2' if _n==2
			replace hi = `m2' + 1.*`se2' if _n==2
			replace lo = `m2' - 1.*`se2' if _n==2				
			twoway (bar e n if n<=2,barwidth(.5)bcol(gs13)ytit(Probability of democratic collapse)saving(h5.gph,replace)) ///
				(rspike hi lo n if n<=2,ylab(0(.02).06)yscale(range(0 .06)) ///
				col(gs1)legend(off)xtit("")tit({bf:Low} system institutionalization)text(.05 1.5 "p<0.01") ///
				xlab(1 "No create party (7)" 2 "Create party (16)")xscale(range(.8 2.2)))
			tab gwf_back create  if ipi>=`c'  ,col
			ttest gwf_back if ipi>=`c',by(create)level(95)
			local  m1=r(mu_1) 
			local se1 = r(sd_1)/(sqrt(r(N_1)))
			local  m2=r(mu_2)
			local se2 = r(sd_2)/(sqrt(r(N_2)))
			replace e=`m1' if _n==1
			replace hi = `m1' + 1.645* `se1' if _n==1
			replace lo = `m1' - 1.645* `se1'  if _n==1
			replace e=`m2' if _n==2
			replace hi = `m2' + 1.645*`se2' if _n==2
			replace lo = `m2' - 1.645*`se2' if _n==2				
			twoway (bar e n if n<=2,barwidth(.5)bcol(gs13)ytit(Probability of democratic collapse)saving(h6.gph,replace)) ///
				(rspike hi lo n if n<=2,ylab(0(.02).06)yscale(range(0 .06)) ///
				col(gs1)legend(off)xtit("")tit({bf:High} system institutionalization)text(.05 1.5 "p<0.05") ///
				xlab(1 "No create party (7)" 2 "Create party (5)")xscale(range(.8 2.2)))
			gr combine h1.gph h2.gph h3.gph h4.gph h5.gph h6.gph, col(2) ysize(7)iscale(.6)
			gr export "$dir\golden\Ch4-Democracy-ttests-by-partytype.pdf",as(pdf)replace
			
			
			local var  = "persparty priormil pres loggdp l1polar ipi"
			foreach v of local var {
				qui reg `v' ld ivdem if year==min
				qui predict r_`v' if e(sample)==1,r
				qui xtset cowcode year
				qui xtsum r_`v' 
				local b = r(sd_b)
				local w = r(sd_w)
				local r = `w'/(`b' +`w')
				di "`v'"
				di `r'
			}
			drop r_* 
			
		**********************
		** Democratic decay **
		**********************
		global d="persparty"
		global ldv = "l1v2x_polyarchy l2v2x_polyarchy l3v2x_polyarchy"
		
		  **** Look at serial correlation ****
		  xtset cowcode year
		  qui xi:xtreg v2x_polyarchy  i.year persparty ld,fe
		  		qui predict e_residuals_1, e
				xtistest e_res, lags(2)
				drop e_res
		  qui xi:xtreg v2x_polyarchy ivdem  i.year persparty ld,fe
		  		qui predict e_residuals_1, e
				xtistest e_res, lags(2)
				drop e_res
		  qui xi:xtreg v2x_polyarchy l1v2x_poly l2v2x_poly  i.year persparty ld,fe
		  		qui predict e_residuals_1, e
				xtistest e_res, lags(2)
				drop e_res
		  qui xi:xtreg v2x_polyarchy l1v2x_poly l2v2x_poly l3v2x_poly i.year persparty ld,fe
		  		qui predict e_residuals_1, e
				xtistest e_res, lags(2)
				drop e_res
	 
				
			replace l1v2x_poly=l1v2x_poly*20
			replace l2v2x_poly=l2v2x_poly*20
			replace l3v2x_poly=l3v2x_poly*20

			* Create party *
			reghdfe v2x_polyarchy create ld ivdem,absorb(cowcode year)cluster(lid)
			est store dem1a
			reghdfe v2x_polyarchy create $ldv ld,absorb(cowcode year)cluster(lid)
			est store dem1b
			* Party personalism *
			reghdfe v2x_polyarchy $d ld ivdem,absorb(cowcode year)cluster(lid)
			est store dem1
			reghdfe v2x_polyarchy $d $ldv ld,absorb(cowcode year)cluster(lid)
			est store dem1c
			
			* Newey, HAC, and PCSE errors *
			qui reghdfe v2x_polyarchy $d $ldv ld,absorb(cowcode year)cluster(lid)
			lincom $d
			qui newey v2x_polyarchy i.cowcode i.year $d $ldv ld,lag(2)force
			lincom $d
			qui ivreghdfe v2x_polyarchy $d $ldv ld,a(cowcode year)bw(3)rob
			lincom $d
			qui xtpcse v2x_polyarchy i.cowcode i.year $d $ldv ld,corr(ar1)het pair
			lincom $d

			* Check IFE *
			regife v2x_polyarchy $d ld ivdem,absorb(cowcode year)ife(cowcode year,1)vce(cluster lid)
			est store dem1d
			regife v2x_polyarchy $d $ldv ld,absorb(cowcode year)ife(cowcode year,1)vce(cluster lid)
			est store dem1e
			regife v2x_polyarchy $d $ldv ld,absorb(cowcode year)ife(cowcode year,2)vce(cluster lid)
			
			label var create "Create party"
			label var l1v2x_poly "Democracy level{sub:t-1}"
			label var l2v2x_poly "Democracy level{sub:t-2}"
			label var l3v2x_poly "Democracy level{sub:t-3}"
			coefplot (dem1a, msymbol(d))(dem1b, msymbol(oh))(dem1, msymbol(t))(dem1c, msymbol(D)) (dem1d, msymbol(Oh)) ///
				(dem1e, msymbol(T)),  order(create persparty ivdem $ldv ld)  ///
				drop(_cons  time) xline(0) msymbol(d) mfcolor(white) grid(glcolor(gs15)) ///
				levels(95 90) legend(lab(3 "Create + initial dem.")lab(6 "Create + lagged dem.") ///
				lab(9 "Personalism + inital dem.")lab(12 "+ Personalism + lagged dem.")  ///
				lab(15 "IFE, Personalism + inital dem.")lab(18 "IFE, + Personalism + lagged dem.") order(3 6 9 12 15 18) ///
				size(small) pos(6) col(2) ring(1)) xsize(2) ysize(2) xlab(-.06(.03).09)  ///
				xtitle("        Coefficient estimate", size(small))  ///
				ciopts(lwidth(thin)) aspectratio(1.1) scale(.75) title(Democratic decay, size(medium) height(2))
			gr export "$dir\golden\T-Decay-Robustness.pdf",as(pdf)replace

			* Check democratic decay against confounders *
			local i =1
			local var = "l2v2xps_party ipi i_pop election i.v2elparlel v2pavote v2paseats ipolar isupdem priormil pres econcrisis gdp lpop oilgas v2pariglef"
			foreach v of local var {
				di "`v'"
				xtset lid year
				qui  reghdfe v2x_polyarchy `v' $ldv $d ld,absorb(cowcode year)cluster(lid)
				lincom $d
				qui nlcom _b[$d],post
				matrix beta =e(b)  
				local b = beta[1,1]
				qui replace beta=`b' if n==`i'
				matrix var = e(V) 
				local se =var[1,1]
				qui replace hi = `b' + sqrt(`se')*1.96 if n==`i'
				qui replace lo = `b' - sqrt(`se')*1.96 if n==`i'
				qui replace hi90 = `b' + sqrt(`se')*1.645 if n==`i'
				qui replace lo90 = `b' - sqrt(`se')*1.645 if n==`i'
				qui replace varname = "`v'" if n==`i'
				local i = `i' +1
			 }
			label define varlab 1  "Party instit." 2 "Initial party inst." 3 "Populism" ///
				4 "Election" 5 "Electoral system" 6 `""Ruling party" "leg. seat share""' ///
				7 `""Ruling party" "vote share""' 8 "Polarization" 9`""Citizen support" "for democracy""' 10 "Prior military" ///
				11 "Presidential" 12 "Economic crisis" 13 "GDP per capita" 14 "Population" 15 "Oil rents" 16 "Right-left ideology",replace
			label values n varlab
			twoway (scatter beta n if n<=16,mcol(blue)yscale(range(-.03 0))yline(-.0187384,lcol(gs4)lpat(dash_dot))) ///
				(rspike hi lo n if n<=16,lw(vthin)lcol(blue)ylab(-0.03(.01)0)) ///
				(rspike hi90 lo90 n if n<=16,lcol(blue)lw(medium)ytitle("{&beta}{sub:Party personalism}", ///
				size(large)height(4))tit(Democratic decay)subtit("Covariate adjustment, dynamic panel model",size(small))   ///
				xtitle(Added covariate,height(-12))yline(0,lpat(dash)lcol(red))xlab(1(1)16,valuelabel angle(90))legend(off))
			gr export "$dir\golden\T-Persparty-DemDecay-covariates.pdf",as(pdf)replace 
			
			* Adjusting for support for democracy does not change estimate for party personalism -- only changing the sample does *
			reghdfe v2x_polyarchy $ldv $d ld if isupdem~=.,absorb(cowcode year)cluster(lid year)
			reghdfe v2x_polyarchy isupdem $ldv $d ld,absorb(cowcode year)cluster(lid year)

		************************
		** Democratic erosion **
		************************
		gen time2 =time^2
			* Check different thresholds *
		forval i = 8(1)15 {
 			di `i'
			local c = `i'/100
			qui reghdfe decline`i' ivdem ld persparty if ld`i'~=1,absorb(year)vce(cluster lid)  /* pooled */
			lincom persparty
			qui reghdfe decline`i' ivdem ld persparty if ld`i'~=1,absorb(cowcode year)vce(cluster lid) /* within */
			lincom persparty
		}	 
		
		 
		* Two-way FE LPM *
		reghdfe decline10 persparty ld if ld10~=1,absorb(cowcode year)cluster(lid)
		tab decline10 if e(sample)==1
		reghdfe decline10 persparty ivdem ld if ld10~=1,absorb(cowcode year)cluster(lid)
 		
		* Check lag  of decline10 the prior year*
		reghdfe decline10 persparty ivdem ld ld10,absorb(cowcode year)cluster(lid)
		tab decline10 if e(sample)==1

		* Check IFE *
		regife decline10 persparty ivdem ld if ld10~=1,absorb(cowcode year)ife(cowcode year,1)vce(cluster lid)
 
		* CRE Probit *
		xi:xthybrid decline10 ld time $d if ld10~=1,cluster(cowcode) vce(cluster cowcode)cre family(binomial)link(probit)se
		tab decline10 if e(sample)==1
		xi:xthybrid decline10 ivdem ld time $d if ld10~=1,cluster(cowcode) vce(cluster cowcode)cre family(binomial)link(probit)se
		est store dem2
		local var =  "ivdem ld time persparty"
		foreach v of local var {
			qui egen xm_`v'=mean(`v') if e(sample)==1,by(cowcode)
		}
		probit decline10 ivdem ld time $d xm_* if ld10~=1,cluster(lid)
		margins,dydx(persparty)
		xtprobit decline10 ivdem ld time $d xm_* if ld10~=1,i(cowcode)vce(cluster cowcode)
		drop xm_*
		
		* Check democratic erosion against confounders *
			local i =1
			local var = "l1v2xps_party ipi v2xpa_popul election i.v2elparlel v2pavote v2paseats ipolar isupdem priormil pres econcrisis gdp lpop oilgas v2pariglef"
			foreach v of local var {
				di "`v'"
				xtset lid year
				qui xi:xthybrid decline10  `v' ivdem ld time $d if ld10~=1,cluster(cowcode) vce(cluster cowcode)cre family(binomial)link(probit)p
				qui lincom  W__persparty
				local b = r(estimate)
				qui replace beta=`b' if n==`i'
				qui lincom  W__persparty
				mat l =r(lb)
				mat h =r(ub)
				local l = l[1,1]
				local h = h[1,1]
				qui replace hi = `h' if n==`i'
				qui replace lo = `l' if n==`i'
				qui replace varname = "`v'" if n==`i'
				local i = `i' +1
			 }
			label define varlab  1 "Party instit." 2 "Initial party inst." 3 "Populism" ///
				4 "Election" 5 "Electoral system" 6 `""Ruling party" "leg. seat share""' ///
				7 `""Ruling party" "vote share""' 8 "Polarization" 9 `""Citizen support" "for democracy""' 10 "Prior military" ///
				11 "Presidential" 12 "Economic crisis" 13 "GDP per capita" 14 "Population" 15 "Oil rents" 16 "Right-left ideology",replace
			label values n varlab
			twoway (scatter beta n if n<=16,mcol(blue)yscale(range(-.01 0.08))yline(1.3029 ,lcol(gs4)lpat(dash_dot))) ///
				(rspike hi lo n if n<=16,lw(vthin)lcol(blue)ylab(0(1)3) ytitle("{&beta}{sub:Party personalism}", ///
				size(large)height(4))tit(Democratic erosion)subtit("Covariate adjustment, initial democracy level panel model",size(small))   ///
				xtitle(Added covariate,height(-12))yline(0,lpat(dash)lcol(red))xlab(1(1)16,valuelabel angle(90))legend(off))
			gr export "$dir\golden\T-Persparty-DemErosion-covariates.pdf",as(pdf)replace 
			drop n
			
			
		* Kernel regression for democratic decline *
		qui reg decline10 persparty ivdem ld time if ld10~=1
		egen cnt =count(year) if e(sample)==1,by(cowcode)
		gen s = e(sample)==1  
		local var = "ivdem ld persparty time"
		foreach v of local var { 
			egen xm_`v' =mean(`v') if s==1,by(cowcode)
		}	
		qui reg decline10 ivdem ld persparty time time2 xm_* if ld10~=1
		lincom persparty
		qui reghdfe decline10 ivdem ld persparty if ld10~=1,a(cowcode year)
		lincom persparty
		 
		krls decline10 ivdem ld persparty time time2 xm_ivdem xm_ld xm_persparty xm_time,d(k)lambda(7.348)
		qui sum persparty
		replace k_pers=k_pers*r(sd)*4  /* probability over four years */
 		ttest k_pers,by(pres)	
		sum k_pers if v2paseats==.
 		sum k_pers if v2paseats!=.
		sum decline10 if s==1  /* 1.4 percent of country-years; so over 4 years: 5.5 percent */

		twoway (lpolyci k_pers v2paseat,bw(15)col(gs12)lpat(dot)) ///
			(lpoly k_pers v2paseat if pres==1,bw(15)lcol(gs6)lpat(solid)xlab(0(20)100)  ///
			legend(lab(1 "All democracies")lab(2 "CI") lab(3 "Presidential") lab(4 "Parliamentary")order(1 3 4)pos(6)ring(1)col(3)) ///
			tit(Party personalism and democratic decline)ytit(Marginal effect of Party personalism)) ///
			(lpoly k_pers v2paseat if pres==0,bw(15) lpat(dash)col(gs1)xtit(Ruling party seat share)  ///
			ytit(Marginal effect of Party personalism)ylab(0(.02).06) )
		gr export "$dir\golden\T-DemDecline-SeatShare.pdf",as(pdf)replace 
		 
		 
		****************************
		** Descriptives by region **
		****************************
		label def region 1 "Eastern Europe and Central Asia" 2 "Latin America and the Caribbean"  ///
			3 "Middle East and North Africa" 4 "Sub-Saharan Africa" 5 "Western Europe and North America" ///
			6 "Asia and Pacific" ,replace
		label val pregion region
		tab pregion if (decline10==1 & ld10~=1) | gwf_back==1
		
		
		*************************
		** Democratic collapse **
		*************************
		drop xm_*
		* Create party *
		reghdfe gwf_back create ld ivdem,absorb(cowcode year)cluster(lid year)
		reghdfe gwf_back create ld $ldv,absorb(cowcode year)cluster(lid year)
		* Party personalism *
		reghdfe gwf_back $d ld ivdem,absorb(cowcode year)cluster(lid year)
		reghdfe gwf_back $d ld $ldv,absorb(cowcode year)cluster(lid year)
		local var ="$d d1 d2 d3 ivdem time"
		foreach v of local var {
			egen xm_`v'=mean(`v') if e(sample)==1,by(cowcode)
		}
		probit gwf_back $d d1 d2 d3 ivdem time xm_*,vce(cluster lid)
		margins,dydx(persparty)
		est store dem3
		drop xm_ivdem
		forval i = 1/3 {
			egen xm_l`i' = mean(l`i'v2x_polyarchy),by(cowcode)
		}
		probit gwf_back $d d1 d2 d3 $ldv time xm_*,vce(cluster lid)
		margins,dydx(persparty)
		* IV-2SLS estimator *
		reghdfe gwf_back v2paind d1 d2 d3 ivdem,absorb(cowcode year)cluster(lid year)
		ivreghdfe gwf_back d1 d2 d3 ivdem (v2paind=create),absorb(cowcode year)cluster(lid year)
		
		*** Marginal effect over time ***
		drop xm_* k_*
		qui reghdfe gwf_back $d time d1 d2 d3 ivdem,absorb(cowcode)cluster(lid)
		lincom $d
		local var = "persparty d1 d2 d3 ivdem time"
		foreach v of local var {
			qui egen xm_`v'=mean(`v'),by(cowcode)
		}
		qui reg gwf_back $d time d1 d2 d3 ivdem xm_*,cluster(lid)
		lincom $d
		qui probit gwf_back $d time d1 d2 d3 ivdem xm_*,cluster(lid)
		margins,dydx($d)	

		krls gwf_back $d d1 d2 d3 ivdem time xm_*,d(k)
		gen kpers=k_pers*.6*3	/* Marginal effect over three years in office of a change from low to high personalism (0.6) */	
		sum kpers
		twoway lpolyci kpers year,ylab(0(.01).04)legend(off)xtit(Year) ///
			ytit("Marginal effect of party personalism" "on democratic collapse",size(small)) ///
			tit(Party personalism increases democratic collapse more in the past decade)
		gr export "$dir\golden\Ch4-Persparty-Collapse-Year.pdf",as(pdf)replace 
			
		* Check democratic collapse against confounders *
		drop xm_*
		gen n=_n
		gen proportional= v2elparlel==1 | v2elparlel==3
		gen mixed =  v2elparlel==2
		qui reghdfe gwf_back $d d1 d2 d3 $ldv,absorb(cowcode year)
		local var ="$d d1 d2 d3 ivdem time"
		foreach v of local var {
			egen xm_`v'=mean(`v') if e(sample)==1,by(cowcode)
		}
		local i =1
		local var="l2v2xps_party ipi i_pop election proportional mixed v2pavote v2paseats ipolar isupdem priormil pres econcrisis gdp lpop oilgas"
		foreach v of local var {
				di "`v'"
				egen mx_`v'=mean(`v'),by(cowcode)
				xtset lid year
				qui probit gwf_back `v' $d d1 d2 d3 ivdem xm_* mx_,cluster(lid)
				lincom $d
				qui nlcom _b[$d],post
				matrix beta =e(b)  
				local b = beta[1,1]
				qui replace beta=`b' if n==`i'
				qui probit gwf_back `v' $d d1 d2 d3 ivdem xm_* mx_,cluster(lid)
				qui lincom $d
				mat l =r(lb)
				mat h =r(ub)
				local l = l[1,1]
				local h = h[1,1]
				qui replace hi = `h' if n==`i'
				qui replace lo = `l' if n==`i'
				qui replace varname = "`v'" if n==`i'
				local i = `i' +1
				drop mx_
			}
			format n %9.0g
			label define varlab 1  "Party instit." 2 "Initial party inst." 3 "Populism" ///
				4 "Election" 5 `""Proportional" "system""' 6 "Mixed system" 7 `""Ruling party" "leg. seat share""' ///
				8 `""Ruling party" "vote share""' 9 "Polarization" 10 `""Citizen support" "for democracy""' 11 "Prior military" ///
				12 "Presidential" 13 "Economic crisis" 14 "GDP per capita" 15 "Population" 16 "Oil rents" ,replace
			label values n varlab,nofix
			twoway (scatter beta n if n<=16,mcol(blue)yscale(range(0 .06))yline( 1.1992,lcol(gs4)lpat(dash_dot))) ///
				(rspike hi lo n if n<=16,lw(vthin)lcol(blue)ylab(-1.5(1)3.5)) ///
				(rspike hi90 lo90 n if n<=16,lcol(blue)lw(medium)ytitle("{&beta}{sub:Party personalism}", ///
				size(large)height(4)) scale(.9)tit(Democratic collapse)subtit("Covariate adjustment, within panel model",size(small))   ///
				xtitle(Added covariate,height(0) )yline(0,lpat(dash)lcol(red))xlab(1(1)16,valuelabel angle(90))legend(off))
			gr export "$dir\golden\T-DemCollapse-covariates.pdf",as(pdf)replace 
	
		* Too much missingness in citizen support from democracy *
		tab gwf_back create if isupdem~=., col
		tab gwf_back create if isupdem==., col 
		reghdfe gwf_back $d ld ivdem time,absorb(cowcode)cluster(lid)
		reghdfe gwf_back $d ld ivdem time if isupdem==.,absorb(cowcode)cluster(lid)
		reghdfe gwf_back $d ld ivdem time if isupdem~=.,absorb(cowcode)cluster(lid)
		reghdfe gwf_back $d ld ivdem time isupdem if isupdem~=.,absorb(cowcode)cluster(lid)
 
		****************************************
		* Democratic decay and erosion results *
		****************************************
		estout dem1 dem2 dem3 using Ch4-Table1.tex,cells(b(star  fmt(%9.3f)) se(par fmt(%9.3f))) ///
			stats(N N_clust) style(tex) replace label starlevels(* 0.10 ** 0.05) title(\label{tab1})
			
		********************
		*** CRE collapse ***
		********************
		drop xm_*
		local var = "persparty i_pop ipolar pres priormil d1 d2 d3 ivdem time time2"
		foreach v of local var {
			egen xm_`v'=mean(`v'),by(cowcode)
		}
		global xm = "xm_persparty xm_d1 xm_d2 xm_d3 xm_ivdem xm_time xm_time2 d1 d2 d3 ivdem time time2"
		qui probit gwf_back $xm persparty,cluster(lid)
		lincom $d
		qui probit gwf_back $xm persparty pres xm_pres,cluster(lid)
		lincom $d
		lincom pres
		qui probit gwf_back $xm persparty priormil xm_priormil,cluster(lid)
		lincom $d
		lincom priormil
		qui probit gwf_back $xm persparty i_pop xm_i_pop,cluster(lid)
		lincom $d
		lincom i_pop
		qui probit gwf_back $xm persparty ipolar xm_ipolar,cluster(lid)
		lincom $d
		lincom ipolar
		qui probit gwf_back $xm persparty pres xm_pres priormil xm_priormi i_pop xm_i_pop,cluster(lid)
		lincom $d
 		lincom pres
		lincom priormil
		lincom i_pop
		sum $d pres priormil i_pop
		margins,dydx($d pres priormil i_pop)
		mat list r(table)
		mat r=r(table)
		forval i=1(1)5 {
			local e = r[1,`i']
			replace e = `e' if n==`i'
			local l = r[5,`i']
			replace lo = `l' if n==`i'
			local h = r[6,`i']
			replace hi = `h' if n==`i'
		}
		twoway (rspike lo hi n if n<=4)  (scatter e n if n<=4,col(gs1)msym(O)yline(0,lcol(red) ///
			lpat(solid)) ytit(Marginal effect)xtit("")legend(off)xscale(range(0.8 4.2)) ///
			xlab(1 "Party personalism" 2 "Presidential" 3 "Prior military" 4 "Populism") ///
			tit(Predictors of democratic collapse))
			
	************************************
	**** By ruling party seat share ****
	************************************
			use pers-use,clear
			gen decline10 = (v2x_polyarchy - ivdem)<=-.1  if  ivdem~=. & v2x_polyarchy~=.
			xtset lid year 
			gen ld10=l.decline10  
			drop b hi* lo*
			gen hi5=.
			gen lo5=.
			gen lo10=.
			gen hi10=.
			global n=3
			
			gen major = v2paseats>=50 if v2paseats~=.
			xi: reghdfe v2x_poly i.major*persparty ld ivdem,a(cowcode year)cluster(lid)
			lincom persparty
			lincom persparty + _ImajXpersp_1
			
			xi: reghdfe v2x_poly i.major*persparty ld ivdem if pres==1,a(cowcode year)cluster(lid)
			lincom persparty
			lincom persparty + _ImajXpersp_1
			xi: reghdfe v2x_poly i.major*persparty ld ivdem if pres==0,a(cowcode year)cluster(lid)
			lincom persparty
			lincom persparty + _ImajXpersp_1
			
			** Plot effects of a move from a low personalism (10pctile) to a high personalism (90pctile) party
			centile persparty if v2paseat~=.,centile(10 90)
		interflex v2x_poly persparty v2paseat ld ivdem,fe(cowcode year)cluster(lid)nbin($n)seed($seed)cutoffs(39.9 49.9)
			mat e=r(estBin)
			gen b=.
			gen se=.
			gen n=_n
			gen xdist=.
			forval i = 1/$n {
  				qui replace xdist=e[`i',1] if n==`i'
  				qui replace b=e[`i',2] if n==`i'
  				qui replace se=e[`i',3] if n==`i'
			}
			replace hi5  = b+1.96*se
			replace hi10 = b+1.65*se
			replace lo5  = b-1.96*se
			replace lo10 = b-1.65*se
			gen rb = round(b,.001)
			graph twoway  (rspike hi5 lo5 n if  n<=$n,lcolor(gs1)lwidth(thin) ///
				xtitle("Ruling party legislative seat share",size(small)height(3)) ///
				yline(0,lp(dash)lcol(red))ytitle("Marginal effect of party personalism", ///
				size(small) height(3))ylab(-.10(.02).00) title("Democratic decay",size(med)))  ///
				(scatter b n if n<=$n,msym(O)mlab(rb) lpattern(solid)lcolor(blue*1.1) ///
				xlab(1 "less than 40%" 2 "40% - 50%" 3 "more than 50%")xscale(range(0.75 $n.25))  ///
				legend(lab(1 "95% CI")lab(2 "Marginal effect")  size(small)pos(7)col(1)ring(0) ///
				order(1 2)))  
			gr export "$dir\golden\Ch4-DemTests-RulingPartySeatShare.pdf",as(pdf)replace 
			drop xdist n b se
		interflex decline10 persparty v2paseat ld ivdem if ld10~=1,fe(cowcode year)cluster(lid)nbin($n)seed($seed)
			mat e=r(estBin)
			gen b=.
			gen se=.
			gen xdist=.
			gen n=_n
			forval i = 1/$n {
  				qui replace xdist=e[`i',1] if n==`i'
  				qui replace b=e[`i',2] if n==`i'
  				qui replace se=e[`i',3] if n==`i'
			}
			replace hi5  = b+1.96*se
			replace hi10 = b+1.65*se
			replace lo5  = b-1.96*se
			replace lo10 = b-1.65*se
			local var ="hi5 hi10 lo5 lo10 b"
			foreach v of local var {
				replace `v'=`v'*.6 /* hi-lo */
			}
			graph twoway  (rspike hi5 lo5 n if  n<=$n,lcolor(gs1)lwidth(thin) ///
				xtitle("Ruling party legislative seat share",size(small)height(3)) ///
				yline(0,lp(dash)lcol(red))ytitle("Marginal effect of party personalism", ///
				size(small) height(0))ylab(-.02(.02).08) title("Democratic declines",size(med)))  ///
				(scatter b n if n<=$n,msym(O) lpattern(solid)lcolor(blue*1.1) ///
				xlab(1 "Lower" 2 "Medium" 3 "Higher")xscale(range(0.75 $n.25))  ///
				legend(lab(1 "95% CI")lab(2 "Marginal effect")  size(small)pos(7)col(1)ring(0) ///
				order(1 2))) (connected  b n if n<=$n,saving(h1.gph,replace)) 
				drop xdist n b se
		interflex gwf_back persparty v2paseat ld ivdem,fe(cowcode year)cluster(lid)nbin($n)seed($seed)
			mat e=r(estBin)
			gen b=.
			gen se=.
			gen xdist=.
			gen n=_n
			forval i = 1/$n {
  				qui replace xdist=e[`i',1] if n==`i'
  				qui replace b=e[`i',2] if n==`i'
  				qui replace se=e[`i',3] if n==`i'
			}
			replace hi5  = b+1.96*se
			replace hi10 = b+1.65*se
			replace lo5  = b-1.96*se
			replace lo10 = b-1.65*se
			local var ="hi5 hi10 lo5 lo10 b"
			foreach v of local var {
				replace `v'=`v'*.6  /* 1 hi-lo */
			}
			graph twoway  (rspike hi5 lo5 n if  n<=$n,lcolor(gs1)lwidth(thin) ///
				xtitle("Ruling party legislative seat share",size(small)height(3)) ///
				yline(0,lp(dash)lcol(red))ytitle("Marginal effect of party personalism", ///
				size(small) height(0))ylab(-.02(.02).08) title("Democratic collpase",size(med)))  ///
				(scatter b n if n<=$n,msym(O) lpattern(solid)lcolor(blue*1.1) ///
				xlab(1 "Lower" 2 "Medium" 3 "Higher")xscale(range(0.75 $n.25))  ///
				legend(lab(1 "95% CI")lab(2 "Marginal effect")  size(small)pos(7)col(1)ring(0) ///
				order(1 2))) (connected  b n if n<=$n,saving(h2.gph,replace))
				drop xdist  b se n
			graph combine h1.gph h2.gph , col(2) xsize(8)	
			gr export "$dir\golden\T-DemTests-RulingPartySeatShare.pdf",as(pdf)replace 

			* Adjust for party personalism *
			qui interflex v2x_poly persparty v2paseat ld ivdem,fe(cowcode year)cluster(lid)nbin($n)
			mat e=r(estBin)
			mat list e
			qui interflex v2x_poly persparty v2paseat v2xpa_popul ld ivdem,fe(cowcode year)cluster(lid)nbin($n)
			mat e=r(estBin)
			mat list e
			qui interflex gwf_back persparty v2paseat ld ivdem,fe(cowcode year)cluster(lid)nbin($n)
			mat e=r(estBin)
			mat list e
			qui interflex gwf_back persparty v2paseat v2xpa_popul ld ivdem,fe(cowcode year)cluster(lid)nbin($n)
			mat e=r(estBin)
			mat list e
			
		*****************************************
		*** Interaction by executive approval ***
		*****************************************
			use "$dir\eap.dta",clear
			drop if year<1990
			egen cid=group(country)
			tsset cid year
			forval i = 1/3 {
				gen l`i'approval = l`i'.approval_not_smoothed
			}
			replace country = "Bulgaria" if country=="Bulgaria_EXEC"
			replace country = "Czech Republic" if country=="Czech Republic_EXEC"
			replace country = "France" if country=="France_EXEC"
			replace country = "Macedonia" if country=="Macedonia_EXEC"
			replace country = "Poland" if country=="Poland_EXEC"
			replace country = "Portugal" if country=="Portugal_EXEC"
			replace country = "Russia" if country=="Russian Federation_EXEC"
			replace country = "Turkey" if country=="Turkey_EXEC"
			replace country = "Ukraine" if country=="Ukraine_EXEC"
			gen cowcode =.
			qui do cowcodes
			sort cowcode year
			merge cowcode year using "$dir\pers-use.dta"
			tab country if _merge==1 & year>1990
			drop if GNB==1
			sort lid year
			gen repeat =year==year[_n-1] if lid==lid[_n-1]
			drop if repeat==1
			drop repeat
			corr l1approval v2paseat v2pavote
			
			gen opolar = l1polar if year==min
			egen ipolar = max(opolar),by(lid)
			
			* Biased sample with nonmissing data on executive approval *
			gen missapproval =  l1approval==. if persparty~=.
			ttest ivdem,by(miss)
			ttest ld,by(miss)
			ttest persparty,by(miss)
			ttest gwf_back,by(miss)

			gen elect_exec = (v2xel_elecparl==1 & pres==0) |  (v2xel_elecpres==1 & pres==1)
			tab elect_exec
 
 			interflex v2x_poly persparty v2paseats ivdem ld elect_exec, ///
				fe(cowcode year leadertime)cluster(lid)nbin(3)cutoffs(39.9 49.9) 
			mat e1=r(estBin)
			interflex v2x_poly persparty v2paseats ivdem ld elect_exec if l1approval~=., ///
				fe(cowcode year leadertime)cluster(lid)nbin(3)cutoffs(39.9 49.9)  
			mat e2=r(estBin)
			interflex v2x_poly persparty l1approval ivdem ld elect_exec if l1approval~=., ///
				fe(cowcode year leadertime)cluster(lid)nbin(3)cutoffs(39.9 49.9)  
			mat e3=r(estBin)
			* Reported models in Table 4.2 *
			forval i=1/3 {
				mat list e`i'
			}
			gen persXseats = persparty*v2paseats
			interflex v2x_poly persparty l1approval ivdem ld elect_exec v2paseat persXseats if l1approval~=.,  ///
				fe(cowcode year leadertime)cluster(lid)nbin(3)cutoffs(39.9 49.9)
				
			* Check whether seat share boost beyond approval is associated with democratic decline *
			reg v2paseats l1approval
			predict res if e(sample),res
			gen devdem = v2x_poly -ivdem
			twoway lpolyci devdem res, legend(off)
			krls devdem res
			reg devdem res,cluster(lid)
			drop res
			reg v2paseats v2pavote
			predict res if e(sample),res
			krls devdem res
			reg devdem res,cluster(lid)
			drop res	
			
			* Check legislative seat share moderation by system type *
			gen lnleader = ln(1+leadertime)
			interflex v2x_poly persparty v2paseats ivdem ld elect_exec lnleader,fe(cowcode year)cluster(lid)nbin(2)cutoffs(51)
			mat list r(estBin)
			interflex v2x_poly persparty v2paseats ivdem ld elect_exec lnleader if pres==0,fe(cowcode year)cluster(lid)nbin(2)cutoffs(51)
			mat list r(estBin)
			interflex v2x_poly persparty v2paseats ivdem ld elect_exec lnleader if pres==1,fe(cowcode year)cluster(lid)nbin(2)cutoffs(51)
			mat list r(estBin)
			interflex v2x_poly persparty v2paseats ivdem ld elect_exec lnleader if v2elparlel==0,fe(cowcode year)cluster(lid)nbin(2)cutoffs(51)
			mat list r(estBin)
			interflex v2x_poly persparty v2paseats ivdem ld elect_exec lnleader if v2elparlel==1 | v2elparlel==3, ///
				fe(cowcode year)cluster(lid)nbin(2)cutoffs(51)
			mat list r(estBin)
			interflex v2x_poly persparty v2paseats ivdem ld elect_exec lnleader if v2elparlel==2,fe(cowcode year)cluster(lid)nbin(2)cutoffs(51)
			mat list r(estBin)
			
			* 5 collapses with nonmissing data on approval *
			tab gwf_back create if l1approval~=.,col
			tab gwf_back create if l1approval==.,col
			interflex gwf_back persparty l1approval ivdem ld,fe(cowcode year)cluster(lid)type(kernel)bw(60)seed($seed)

		*******************************
		**** Coups and power-grabs ****
		*******************************
		use pers-use,clear
		drop if persparty==.
		* Update for collapse in 2020 *
		recode gwf_back (0=1) if year==2020 & (country=="Mali")
		* GNB 2012: non-personalist leader dies a naturual death 9 Jan 2012; 
		* interim president takes over and their is a coup in April 
		* drop the interim leader who is unelected *
		drop if GNB==1
 		gen backtype=""
		replace backtype ="rebellion" if country=="Afghanistan" & year==2019
		replace backtype ="powergrab" if country=="Armenia" & year==1994
		replace backtype ="rebellion" if country=="Azerbaijan" & year==1993
		replace backtype ="powergrab" if country=="Bangladesh" & year==2014
		replace backtype ="powergrab" if country=="Benin" & year==2019
		replace backtype ="coup" if country=="Bolivia" & year==2019
		replace backtype ="coup" if country=="Burundi" & year==1996
		replace backtype ="powergrab" if country=="Burundi" & year==2010
		replace backtype ="rebellion" if country=="Central African Republic" & year==2003
		replace backtype ="coup" if country=="Egypt" & year==2013
		replace backtype ="powergrab" if country=="Guinea Bissau" & year==2002
		replace backtype ="powergrab" if country=="Haiti" & year==1999
		replace backtype ="powergrab" if country=="Hungary" & year==2018
		replace backtype ="coup" if country=="Madagascar" & year==2009
		replace backtype ="coup" if country=="Mali" & year==2012
		replace backtype ="coup" if country=="Mali" & year==2020
		replace backtype ="coup" if country=="Mauritania" & year==2008
		*replace backtype ="coup" if country=="Nepal" & year==2002
		replace backtype ="powergrab" if country=="Nicaragua" & year==2016
		replace backtype ="coup" if country=="Niger" & year==1996
		replace backtype ="powergrab" if country=="Niger" & year==2009
		replace backtype ="coup" if country=="Pakistan" & year==1999
		replace backtype ="powergrab" if country=="Peru" & year==1992
		replace backtype ="rebellion" if country=="Republic of Congo" & year==1997
		replace backtype ="powergrab" if country=="Russia" & year==1993
		replace backtype ="powergrab" if country=="Serbia" & year==2018
		replace backtype ="rebellion" if country=="Sierra Leone" & year==1997
		replace backtype ="powergrab" if country=="Sri Lanka" & year==2010
		replace backtype ="powergrab" if country=="Sri Lanka" & year==2020
		replace backtype ="coup" if country=="Thailand" & year==1991
		replace backtype ="coup" if country=="Thailand" & year==2006
		replace backtype ="coup" if country=="Thailand" & year==2014
		replace backtype ="powergrab" if country=="Turkey" & year==2016
		replace backtype ="powergrab" if country=="Ukraine" & year==2012
		replace backtype ="powergrab" if country=="Venezuela" & year==2005
		replace backtype ="powergrab" if country=="Zambia" & year==1996
		gen coupback = backtype=="coup"  if gwf_back~=.
		gen grabback = backtype=="powergrab"  if gwf_back~=.	
		tab coupback grabback if gwf_back==1,row col
 			gen e=.	
			gen hi=.
			gen lo=.
			gen n =_n
 			gen hipers = persparty>.545 if persparty~=.			
			tab gwf_back hipers,col
			ttest gwf_back,by(hipers)
			local  m1=r(mu_1) 
			local se1 = r(sd_1)/(sqrt(r(N_1)))
			local  m2=r(mu_2)
			local se2 = r(sd_2)/(sqrt(r(N_2)))
			replace e=`m1' if _n==1
			replace hi = `m1' + 1.96* `se1' if _n==1
			replace lo = `m1' - 1.96* `se1'  if _n==1
			replace e=`m2' if _n==2
			replace hi = `m2' + 1.96*`se2' if _n==2
			replace lo = `m2' - 1.96*`se2' if _n==2				
			twoway (bar e n if n<=2,barwidth(.5)bcol(gs13)ytit(Probability of democratic collapse)saving(h1.gph,replace)) ///
				(rspike hi lo n if n<=2,ylab(0(.01).04)col(gs1)legend(off)xtit("")tit(Democratic collapse) ///
				xlab(1 "Low (6)" 2 "High (29)")xscale(range(.8 2.2)))
			tab coupback hipers,col
			ttest coupback,by(hipers)
			local  m1=r(mu_1) 
			local se1 = r(sd_1)/(sqrt(r(N_1)))
			local  m2=r(mu_2)
			local se2 = r(sd_2)/(sqrt(r(N_2)))
			replace e=`m1' if _n==1
			replace hi = `m1' + 1.96* `se1' if _n==1
			replace lo = `m1' - 1.96* `se1'  if _n==1
			replace e=`m2' if _n==2
			replace hi = `m2' + 1.96*`se2' if _n==2
			replace lo = `m2' - 1.96*`se2' if _n==2				
			twoway (bar e n if n<=2,barwidth(.5)bcol(gs13)ytit("Probability of coup collapse")saving(h2.gph,replace)) ///
				(rspike hi lo n if n<=2,ylab(0(.005).02)col(gs1)legend(off)xtit("")tit(Military coup collapse) ///
				xlab(1 "Low (1)" 2 "High(11)")xscale(range(.8 2.2)))
			tab grabback hipers,col
			ttest grabback,by(hipers)
			local  m1=r(mu_1) 
			local se1 = r(sd_1)/(sqrt(r(N_1)))
			local  m2=r(mu_2)
			local se2 = r(sd_2)/(sqrt(r(N_2)))
			replace e=`m1' if _n==1
			replace hi = `m1' + 1.95* `se1' if _n==1
			replace lo = `m1' - 1.95* `se1'  if _n==1
			replace e=`m2' if _n==2
			replace hi = `m2' + 1.95*`se2' if _n==2
			replace lo = `m2' - 1.95*`se2' if _n==2				
			twoway (bar e n if n<=2,barwidth(.5)bcol(gs13)ytit(Probability of power-grab collapse)saving(h3.gph,replace)) ///
				(rspike hi lo n if n<=2,ylab(0(.005).02)col(gs1)legend(off)xtit("")tit(Power-grab collapse) ///
				xlab(1 "Low (3)" 2 "High (15)")xscale(range(.8 2.2)))
			gr combine h1.gph h2.gph h3.gph, xsize(4)ysize(2)col(3)subtit(Party personalism,pos(6))
			gr export "$dir\golden\Ch4-Powergrab-Coup-collapse-ttests.pdf",as(pdf)replace 
			
			*** Democratic collapse via rebellion ***
			list country current_leader year backtype hipers persparty if gwf_back==1 & coupback==0 & grabback==0
			** 1 personalist party score just below the median but other 4 are in the 90th percentile or above **
			
			*** Change in democracy prior to collapse events ***
			gen rebelback = backtype=="rebellion"
			drop max*
			egen maxyr = max(year),by(lid)
			gen l1devdem = l1v2x_poly-ivdem
			gen maxdev = l1devdem if year==maxyr
			
			ttest maxdev if rebelback==0 & coupback==0,by(grabback)
			ttest maxdev if rebelback==0 & grabback==0,by(coupback)
			ttest maxdev if grabback==0  & coupback==0,by(rebelback)

			*** 
			gen time = year-1990
			gen xtime = time/10
			gen xivdem=ivdem*10
			local var = "ld xivdem persparty xtime"
			foreach v of local var {
				egen xm_`v'=mean(`v'),by(cowcode)
			}
			xthybrid grabback xivdem ld persparty xtime,cluster(cowcode)vce(cluster cowcode)family(binomial)link(probit)p
			est store coup1
			probit grabback ld xivdem persparty xtime xm_*,cluster(lid)
			est store coup2
			margins,dydx(persparty)
			xthybrid coupback xivdem ld persparty xtime,cluster(cowcode)vce(cluster cowcode)family(binomial)link(probit)p
			est store coup3
			probit coupback ld xivdem persparty xtime xm_*,cluster(lid)
			est store coup4
			margins,dydx(persparty)
			
			coefplot (coup1, msymbol(D)mfcolor(white))(coup3, msymbol(P)mfcolor(gs1)), ///
				keep(W__persparty W__xivdem W__ld W__xtime )  order(W__persparty W__xivdem W__ld W__xtime ) ///
				drop(_cons  B__*) xline(0)  grid(glcolor(gs15)) levels(95 90) ///
				coeflabels(W__persparty  = `""Ruling    " "party     " "{bf:personalism}""' ///
				W__xivdem= `""Initial    " "democracy" "level     ""' W__ld = "Democracy age" ///
				W__xtime ="Time trend") legend(lab(3 "Power-grab collapse")lab(6 "Coup collapse.")  ///
				order(3 6)size(small) pos(6) col(2) ring(1)) xsize(2) ysize(2) xlab(-1(1)3) ///
				xtitle("        Coefficient estimate", size(small))ciopts(lwidth(thin)) aspectratio(1.1) ///
				scale(.95) title(Power-grab and coup collapses, size(medium) height(2)) ///
				note("Correlated random effects estimator:" "within-estimates reported; between estimates omitted.",size(vsmall)pos(6))
	 		gr export "$dir\golden\T-Powergrab-Coup-collapse-estimates.pdf",as(pdf)replace 

			
	erase .pdf
	forval i=1(1)6 {
		qui erase h`i'.gph
	}
			
	****************** The End *****************
			
	log close

	 
	
 	
